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A longstanding goal of nonequilibrium statistical mechanics has been to extend the conceptual 
power of the Boltzmann distribution to driven systems. We report some new progress towards 
this goal. Instead of writing the nonequilibrium steady-state distribution in terms of perturbations 
around thermal equilibrium, we start from the linearized driven dynamics of observables about 
their stable fixed point, and expand in the strength of the nonlinearities encountered during typical 
fluctuations away from the fixed point. The first terms in this expansion retain the simplicity of 
known expansions about equilibrium, but can correctly describe the statistics of a certain class of 
systems even under strong driving. We illustrate this approach by comparison with a numerical 
simulation of a sheared Brownian colloid, where we find that the first two terms in our expansion 
are sufficient to account for the shear thinning behavior at high shear rates. 


I. INTRODUCTION 

For over a century, the formalism of equilibrium statis¬ 
tical mechanics has provided a powerful means to explain 
how the macroscopic properties of many-body systems at 
thermal equilibrium arise from the microscopic interac¬ 
tions that occur among their constituent parts. The cen¬ 
terpiece in this approach is the Boltzmann distribution, 
which posits that the probability of observing an equili¬ 
brated system with energy e in microstate x at tempera¬ 
ture ksT = 1/(3 is proportional to the so-called “Boltz¬ 
mann weight” pbz(x) oc exp[— (3e{x)\. The key assump¬ 
tion used in deriving the Boltzmann distribution is that 
the system has spent an “ergodically” long time in con¬ 
tact with its surrounding heat bath, so that the combined 
set-up of bath and system is equally likely to be in any ar¬ 
rangement that is allowed by conservation of energy. As a 
result, a quantity evaluated for the system at one instant 
in time (namely, e(x)) can immediately be translated into 
a probability of occurrence for the state x. This micro¬ 
scopic result can be coarse-grained to yield the probabil¬ 
ity of observing the system in a given macroscopic state 
X , defined as a set of microstates that share the same 
values of some observable properties. The coarse-grained 
probability pbz(A') oc can then be written in 

terms of a free energy F(X) = —fc s Tln e _ ^ e ( x )). 

Once time-varying fields drive the system from equilib¬ 
rium by changing the energies e(x, t ) on timescales com¬ 
parable to the system’s relaxation time, the story must 
necessarily become more complicated. In the arbitrary 
nonequilibrium scenario, the probability of being at a 
given location in phase space at time t clearly can de¬ 
pend strongly on where the system was at some earlier 
moment. There is, however, a tempting special case to 
consider even when the Boltzmann distribution does not 
apply: in circumstances where e(x , t) is periodic, the sys¬ 
tem may still ergodically lose its memory of initial con¬ 
ditions after enough time in contact with the fluctuating 
bath. In such a case, it is reasonable to consider whether 
the Boltzmann distribution admits a generalization, in 
which the probability of observing the system in a par¬ 
ticular state after the memory of the initial state is lost 


can still be related exactly to some function of thermo¬ 
dynamic observables. 

Yamada and Kawasaki answered this question in the 
affirmative almost fifty years ago, when they derived 
an effective partition function for a generic nonequilib¬ 
rium steady state in terms of correlations in the currents 
of conserved quantities passing through the system [T], 
launching a fruitful field of research on the exact micro¬ 
scopic distribution in driven steady states 00. It is now 
well known that the simplicity of the Boltzmann weight 
cannot be reproduced in the microscopic probability dis¬ 
tribution for an arbitrary driven steady state, which de¬ 
pends in general on all orders of the time correlations in 
the currents over the system’s past history. 

Thus any attempt to uncover simple principles beneath 
the statistics of nonequilibrium steady states must begin 
by specifying a particular regime of applicability, where 
certain simplifying approximations become valid. The 
near-equilibrium regime was the first to receive careful 
study, leading to an elegant representation of the steady- 
state distribution in terms of the dissipation due to exter¬ 
nally imposed thermal gradients, chemical potential gra¬ 
dients, and velocity fields [Bj. This “McLennan ensem¬ 
ble” applies to a wide range of near-equilibrium systems, 
and can be obtained through a variety of independent 
routes (cf. [7]). Most recently, it has been shown that 
this form can be derived in an especially transparent way 
from the assumption of “microscopic reversibility,” which 
holds for a wide class of physical systems 0119]. 

When external drives become arbitrarily strong, the 
steady-state distribution of observables in a generic phys¬ 
ical system no longer follows this form. But some intution 
about this regime can be built up around a special case, 
where a simple expression for the distribution does exist 
for arbitrarily strong driving, with the general expression 
written in terms of perturbations about this case. In this 
paper, we derive such an expansion about the case where 
fluctuations in the instantaneous dissipative current in 
the driven system obey a linear overdamped Langevin 
equation with additive white noise. In this scenario, a 
form essentially equivalent to the McLennan ensemble 
can be shown to hold regardless of the drive strength, 
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and to be compatible with large departures from the pre¬ 
dictions of linear-response theory. 

This approach bears some resemblances to macroscopic 
fluctuation theory, which derives various statistical prop¬ 
erties of far-from-equilibrium macroscopic systems from 
some minimal restrictions on the form of the macroscopic 
dynamics and the character of the current fluctuations 
jlOj . In particular, a simple form for the driven steady- 
state distribution is obtained that looks very similar to 
the McLennan ensemble, but holds for arbitrarily strong 
driving with arbitrary nonlinearities m- This result is 
obtained through a decomposition of the dissipative cur¬ 
rent into “symmetric” and “asymmetric” parts, which 
purchases a broader range of applicability by sacrific¬ 
ing the immediate physical meaning of the terms in the 
original McLennan form. We do not take advantage of 
this decomposition, but instead seek a representation of 
the distribution entirely in terms of the equilibrium free 
energy and the statistics of the bare externally applied 
work. 

We start in section [Tl] by deriving an exact expres¬ 
sion for the steady-state distribution of an arbitrary ob¬ 
servable in a system with microscopic reversibility whose 
distribution relaxes exponentially to a unique stationary 
state. Then in section |III| we evaluate this expression 
in our special case, and compute the first correction in 
the perturbation expansion. In sections |IV| and |IVD| we 
apply this expansion to a sheared colloid, as a specific 
example of a strongly driven system. Section [V] gives a 
quantitative comparison with a numerical simulation of 
the sheared colloid, showing that the McLennan-like part 
is sufficient to reproduce the qualitative non-linear re¬ 
sponse behavior, and that the first correction term gives 
good quantitative agreement deep into the shear thinning 
regime. Finally, in section VI we describe the thermody¬ 
namic intuition that can be extracted from this new form 
for the distribution, and define some avenues for further 
investigation. 


II. DERIVATION OF DISTRIBUTION 
A. Microscopic Reversibility 

Consider a generic physical system coupled to a large 
“heat bath” of inverse temperature (3 = 1/ksT , but oth¬ 
erwise isolated, so that the total energy of the system 
plus bath can only be changed by external manipulation 
of a specified set of “control parameters” A. These con¬ 
trol parameters directly affect only the system energy 
e(x , A) and not the bath. As usually assumed in statisti¬ 
cal mechanics, the potential energy of the interactions of 
system components with bath components is taken to be 
small compared with the energy in the system, so that 
the total energy can be cleanly divided between the en¬ 
ergy “in the system” and the energy “in the bath.” The 
bath volume is also held fixed, and the system volume is 
only allowed to vary if it is chosen as one of the control 


parameters A. 

The microstate x of the system evolves according to a 
stochastic process, due to its interactions with the fluctu¬ 
ating heat bath. If the combined setup is modeled with 
classical mechanics, x specifies generalized positions and 
momenta of all degrees of freedom in the system (not 
including the bath degrees of freedom), x can also be 
taken to be a discrete microstate label, with jumps be¬ 
tween microstates governed by a rate matrix W xx > (A(t)). 
In this case, thermodynamic consistency is imposed by 
requiring that the rate matrix should eventually equili¬ 
brate the system to the Boltzmann distribution if the A’s 
are held fixed. 

We keep the system out of equilibrium by varying the 
A’s according to a given protocol, which gives the mi¬ 
crostate energies e(x , A(£)) an explicit time dependence. 
When we average over trajectories below, we will always 
assume that the A(f) protocol is held fixed. 

Building on the work of C. Jarzynski |T2], G. Crooks 
has shown that a number of important results concerning 
the nonequilibrium behavior of such a system can be de¬ 
rived from what he calls the “microscopic reversibility” 
condition [5[: 


p R [x*(At-t)\x* 2 ] = e -pQ F[x (t)} 
p F [x(t)\xi] 


(1) 


Here x(t) is a system trajectory of duration At and the 
heat Q F [x(t)] is the energy transferred from the system to 
the bath over the course of that trajectory. The left-hand 
side contains the probability of taking the time-reversed 
path x*(At — t) given a starting state x%, divided by 
the probability of taking the forward path x(t) given the 
starting state x\. The * indicates time-reversal of the 
microstate (changing the signs of all the momenta in the 
classical model). The R and F subscripts refer to the 
driving protocol A(£): the F probabilities and heat are 
computed with the protocol forward from time — At to 
time 0, and the R probabilities have it running in reverse 
from time 0 to time —At. In the sheared colloid we will 
analyze below, the R and F quantities are computed with 
the shear applied in opposite directions. 

Crooks showed that condition ([I]) holds for stochas¬ 
tic dynamics of discrete states under the thermodynamic 
consistency requirement given above [5]. This relation 
can also be derived directly from the time-reversibility 
and phase space conservation of Hamiltonian dynamics, 
if the combined system-plus-bath setup can be treated as 
a closed Hamiltonian system driven by an explicit time- 
dependence in the system Hamiltonian (cf. [T3] for the 
basic approach, although a slightly different result is dis¬ 
cussed there). The expression can be generalized to allow 
particle fluxes into and out of chemical baths, in which 
case an extra term involving chemical potentials must be 
added to the Qf in the exponent mm- 
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B. Coarse Graining 

Now we group the system microstates x according to 
some observable properties, following the approach de¬ 
tailed in m- We can give each group a unique name, 
which will be generically represented by the capital letter 
X. Below, in the analysis of colloidal steady states, X 
will stand for a single number (the mean shear stress at 
the wall) characterizing the configuration of all the col¬ 
loidal particles. In general, A' represents a label for some 
group of microstates, which could be picked out using any 
definite procedure. X* will refer to the group consisting 
of the time-reversed versions x* of all the microstates in 
X. 

The coarse-grained version of equation 0 involves 
transition probabilities among the different groups la¬ 
beled by different X values. This probability will in gen¬ 
eral depend on how the initial state was prepared, since 
different protocols will give rise to different probability 
distributions of microstates within the macrostate. In 
this paper, however, we consider only systems with finite 
relaxation times, and our goal is to analyze the steady 
state that is reached after all correlations between cur¬ 
rent X values and the initial conditions have died out. In 
such a case, the choice of initial distribution becomes ir¬ 
relevant to the steady-state statistics, and we can choose 
the distribution that gives the most useful form for the 
final steady-state distribution of X. 

We will call the initial distributions for the for¬ 
ward and reverse trajectories pi(x) and p 2 (x) respec¬ 
tively, and choose them to be Boltzmann distributions 
p(x) = exp[—/3(e(x, A(t)) — F(X, A(f)))] over the mi¬ 
crostates in the respective macrostates A’i and A|. 
F(X, X(t)) = -k B TlnJ2 xe x exp[-/3e(a;, A(t))] gives the 
proper normalization for a distribution defined only over 
microstates x in macrostate X (where the sum becomes 
an integral in the classical case). To investigate steady- 
state behavior, we must vary these fields periodically, 
and consider trajectories whose duration At is an integer 
multiple of the period. This implies that e(x, A(—At)) = 
e(x, A(0)) and F(x,X(—At)) = F(x, A(0)), so for both 
cases we can drop the time-dependence from the nota¬ 
tion. 

With these distributions in hand, we multiply equa¬ 
tion (fTI) by the denominator of the left-hand side and by 
P2(x 2b then integrate over all trajectories x(t) connect¬ 
ing states X\ in a given macrostate X\ to states x 2 in 
another macrostate A2: 

[ V[x(t)\p 2 (x* 2 )p R [x*(At-t.)\x* 2 \ 

Jx i-*-X 2 

= f V[x(t)]e^ l3QF ^ x ^ t) ^ pi(x{)p F [x(t)\xi\. 

Jxi^x 2 Pl{Xl) 

( 2 ) 

We can simplify this expression by introducing the 


macroscopic transition probabilities 

ttf(AT ->• X 2 ) = [ V[x{t)]pi{xi)p F {x(t)\xi] ( 3 ) 

kr(X 2 ->Xf)= f V[x(t)]p2(x%)p R [x*(At - t)\x 2 ] 

JX 2 ^fX 2 

(4) 

which are defined as the sums of the probabilities of 
all microtrajectories x(t) that accomplish the indicated 
macroscopic transition, given that the system begins in 
the indicated probability distribution over microstates 
(x\ or x 2 ) in the starting macrostate (A'i or Af). Us¬ 
ing these definitions to normalize the distributions over 
trajectories in the integrands in equation ([ 2 ]), we find: 

«n{X 2 -+ Xl ) 

x 7 Tf(Ai —> A 2 ). (5) 

The average {•)x 1 -*x 2 is over a ll trajectories x(t) connect¬ 
ing some microstate x± in X\ (chosen from distribution 
Pi) to some other microstate X2 in A' 2 . 

Now we insert the explicit expressions for pi and p 2; 
to find 

7Tj?.(A| -> X*) P[F(X 2 )-F(X 1 )] 

ttf(Ai —y X 2 ) 

X / e -0KM*(*)]-<*l)+e(x2)]\ _^ x 

^ 2 ( 6 ) 

We have dropped the * in the arguments of e and F, 
because the energy is symmetric under reversal of the 
signs of the momentum coordinates. Now we note that 
by conservation of energy, the work done on the system 
by the variation of the control parameters A(f) over a 
trajectory x(t) from x\ to X 2 satisfies Wf = QF[x{t)\ + 
e(x 2 ) — e(xi). We can thus rewrite the above expression 
as 

Tr(A 2 * ->■ A£) r[f{x 2 )-f(x 1 )) 

ff{Xi A 2 ) 

x {e -fWrW)] )x ^ Xa . ( 7 ) 

The motivation for our choice of the Boltzmann distribu¬ 
tion for pi and P2 is precisely because it replaces the heat 
in the exponent with the work. Since the work is zero 
in the undriven case, where the control parameters are 
fixed, this choice splits the right-hand side cleanly into 
two factors, an equilibrium contribution and a nonequi¬ 
librium correction. Other choices could be made for these 
distributions, and they would generate valid alternative 
forms of the steady-state distribution. 

We can use this expression to compare the probabilities 
of forward transitions from one state Xq to two different 
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states X\ and X 2 


In np (Xo -^l) 

Tf(X 0 —> X 2 ) 


=P[F{X 2 ) - F(X 1 )} 


, (e-^^xo^x, 

{ e - 0 w F Mt )]) Xo ^ X2 

k *«(*?-> *5) 
f r (X* ->• Xq) 


( 8 ) 


This expression is interesting in its own right, showing 
how the relative probabilities of different future possi¬ 
bilities depend not only on the free energies of the pos¬ 
sible future states, but also on the work done on the 
way there and on the “durability” measure contained in 
the reverse probabilities (see [18] for a detailed analysis 
of the physical implications of the corresponding terms 
in a closely related expression). We now specialize to 
a system in which temporal correlations decay exponen¬ 
tially (or faster) with a finite relaxation time r. In this 
case, 7r (Xj —► Xf) must become independent of X, for 
At r, and simply be equal to the steady-state proba¬ 
bility p ss (Xf) of being in the final state: 


ln^M=Pi F ( X 2)- F ( x i)} 

Pss( A 2] 

<e-^[W)]) M2 


+ In 


p£(Xq) 

p£(*5)' 


(9) 


From this equation we can directly extract p p (X) up to 
an overall constant A f that is independent of X : 


lnp ss (X) = - pF[X) - lim ln(e~^ w ) Xo ^ x + AT 

At—too 

( 10 ) 


where we have dropped the F’s from p ss and W because 
we will only be considering “forward” quantities from 
now on. 

Even though the work W in the exponent increases 
without bound as At —» oo, this expression is well-defined 
if one takes care to perform the average before taking the 
limit. For exponentially relaxing systems, {e~P w ) Xo -> x 
converges to a finite value as At increases beyond the 
relaxation time r. This can be seen by considering the 
behavior of the LHS of equation (J7| under these condi¬ 
tions: 7Tp(Xi —► X 2 ) and 7 Tr(X| —> X*) approach the 
finite limiting values ofp^(X 2 ) andp^(X*), respectively, 
for At>r. 

The individual cumulants of W do in general be¬ 
come infinite as At —> oo, however. To obtain a 
useful series representation of the exponential aver¬ 
age term, we can add the constant, finite quantity 
limAt-Kx> \a.{e~ 0w ) Xo —r x 0 to the RHS, and compensate 
by adjusting the normalization constant A/”. Represent¬ 
ing both exponential average terms by their cumulant 
expansions, we obtain: 


lnp„(X) = - 0F(X) - (£) 

' n =1 

(id 

n= 1 ' 

We can now rewrite this expression in terms of the 
finite differences between the cumulants for the trajec¬ 
tories ending in X and the corresponding cumulants for 
the trajectories ending in Xo: 

A(ITT(X) = Km ((W n ) c Xo ^ x - <tn W„) 

At—too 

= lim ((r)^-(r)Uo)- (12) 

Since we have assumed that the system’s state becomes 
decorrelated from its past history after a finite time t, 
the expression in parentheses becomes independent of At 
and remains finite as we take At —^ oo. In the second 
line, the ss —► X averages are over the trajectories whose 
initial conditions are sampled form the steady state dis¬ 
tribution, and that end in state X. The equality follows 
from the fact that the contribution of the initial relax¬ 
ation of the system from X 0 to the steady state is the 
same for both terms on the RHS of the first line. 

Thus we obtain: 

In Pss (X) = - (3F(X) - Y, ( -^ L A(W n ) c (X) +J\f. 

n\ 

71=1 

(13) 


Note that if X is an observable quantified by a continuous 
parameter, p ss (X) can be regarded as a probability den¬ 
sity function for X. In this case, the LHS should strictly 
be written as ln[p ss (X)<5X], where a microstate counts as 
part of macrostate X if the value of the observable lies 
within SX of X'. But the <5X can be chosen to be the 
same for all X, and rolled into the overall constant A f. 

Aside from the coarse-graining step, the general pro¬ 
cedure we have followed thus far for extracting a steady- 
state distribution from microscopic reversibility and ex¬ 
panding in cumulants resembles the derivation by Ko¬ 
matsu et. al of an expansion about equilibrium in the 
strength of the driving field [8]. Expressions like (13) ob¬ 
tained in this way contain cumulants of all orders, and 
only provide new physical insight if the series converges 
rapidly. Komatsu et al. show how a different way of writ¬ 
ing the microscopic reversibility assumption leads to an 
expansion about equilibrium that converges particularly 
rapidly, and is accurate to second order in the strength 
of the drive without including any cumulants of higher 
order than 1 (8]. We are taking a different approach, 
treating equation (131 as an expansion in the size of the 
nonlinearities in the coarse-grained equations of motion 
near the steady state. The coarse-graining allows the pa¬ 
rameters of the linearized dynamics to have a non-trivial 
dependence on the strength of the drive, so that the pre¬ 
dictions of near-equilibrium linear response theory can 
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break down while our expansion parameter is still near 
zero. In the next section, we lay out the details of our 
proposed expansion in terms of a Langevin model for the 
coarse-grained dynamics. Then we will use a simulated 
sheared colloid to illustrate a concrete case where the new 
expansion converges rapidly in a strongly driven system. 

Writing the distribution in terms of the A(W") c ’s is 
convenient for the initial presentation of the theory, but 
for comparison with the colloid simulation, we need to 
work with finite Af’s. Thus we define 

(w n ) c At (x) = (w n y ss ^x (i4) 

«A(r) c (i) + (r)^ 0 

where the averages are all over trajectories of length At, 
and the approximation holds for At r. Since the sec¬ 
ond term in the second line is independent of X , replac¬ 
ing A(IT") C (A) with {W n ) c At (X) in any expression for 
the steady-state distribution only affects the normaliza¬ 
tion. In terms of this new quantity, we thus obtain the 
alternative form 

In Pss(X) « - pF(X) - J2 ^-{W n ) c At {X) +J\f. 

Z —' n! 

n—1 

(15) 


III. PERTURBATION ANALYSIS 

Non-equilibrium steady states are distinguished from 
equilibrium states by the existence of non-zero mean cur¬ 
rents J ss of some conserved quantities (mass density, mo¬ 
mentum, charge, etc.), so that the external field E conju¬ 
gate to J does work on a system of volume V at a mean 
rate W = VEJ SS . In this section, we show that the cumu- 
lant differences A (W n ) c for n > 1 all vanish in the spe¬ 
cial case where a set of coarse-grained variables including 
the relevant J's can be found whose combined dynamics 
are described by a linear overdamped Langevin equation 
with additive white noise. This assumption of linear¬ 
ity does not imply a restriction to the linear response 
regime of near-equilibrium thermodynamics, as we show 
in the following subsection by computing the deviation of 
J SS (E) from its linear-response form in terms of the pa¬ 
rameters of the linear model. In the final subsection we 
will examine small perturbations around this regime to 
see how the cumulant differences grow as nonlinearities 
are introduced. 

For some of the calculations, we will assume that the 
observable X whose probability is being computed is one 
of the currents J. This assumption can be relaxed by a 
simple change of variables, whose Jacobian will be ab¬ 
sorbed into the equilibrium term F(X), as long as the 
dynamics of the current still satisfy the given require¬ 
ments. 


A. Linear Regime 

To compute the conditional averages (•) ss _ i .x, we start 
by writing down an equation of motion for the observ¬ 
ables in our system, such that the trajectories of the ob¬ 
servables for a given realization of the noise can be found 
by solving the equation and applying the final condition 
that the trajectory ends in X at t = 0. We will allow for 
an arbitrary number of observables, contained in a vector 
X, but will allow for only one current J with steady-state 
value J ss that is responsible for the steady-state work. 
This restriction can be relaxed without affecting the fi¬ 
nal result, but it simplifies the intermediate notation. 

We start with the case of a linear Langevin equation 
for the fluctuation dynamics: 

X = A(E)X + B(E)£(t) (16) 

where we have defined X such that X = 0 is the most 
probable value in the steady state. We will choose the 
first element of X to contain the current, so that X\ = 
J — J ss . The matrices A and B are constant in time, 
but depend on the strength of the external field E. £(i) 
is a vector of Gaussian white noise, characterized by its 
mean (£*(£)) = 0 and two-point function = 

S(t - t')Sij. 

With the ansatz X = e At [X(0) + f(t)], we find 

AX + e At J = AX + B^(t) 

d i=e~ M m) 

r° 

f (t) = - J dt'e~ At B£(t') (17) 

so that X = e j4 *X(0) - f° dt’ e^*"*'* B^(t') (keeping in 
mind that t < 0, since we are specifying the final condi¬ 
tion X(0)). 

Since the first observable X\ is equal to the deviation 
J — J ss from the mean steady state current, the work 
done over a given trajectory X(t) is 

W = VE [ dt[X 1 {t) + J ss ] (18) 

J — At 

= VE f dt [(e A %X,(0) + (e At )ijfj(t) + J ss ] , 

(19) 

where we are implicitly summing over all the j’s, using 
the Einstein summation convention. 

We can now compute A(IT)(A,;(0)) as a function of 
any one of the parameters X;(0), by averaging W over 
paths that start in the steady state at time t = — At —> 
—00, and end at specified values of Aj(0) at time t = 0: 

A(W)(Aj(0)) = <W} ss ^ Xi(0) - (W) ss ^o 

= ve f dtie^yyxy o)) Xi(0) . 

J —00 

(20) 
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To obtain the higher cumulants, we first use the fact 
that the solution X(t) obtained above is a sum of inde¬ 
pendent Gaussian random variables, to conclude that the 
steady-state distribution for this linear relaxation is itself 
a Gaussian, which can be written as 


p ss (X) (xexpj-^ C^XjXk 


( 21 ) 


where C is the covariance matrix. 

We can similarly show that the work distribution is 
Gaussian. Since all cumulants beyond the second are 
zero for a Gaussian distribution, we need only compute 
the contribution from the variance: 

A {W 2 Y =(W /2 )s S ->.Y, - (W /2 )ss-> 0 

=V 2 E 2 f dt' f dt (e At )ij(e At )ik 

J — oo J — oo 

x «X,(0)X fc (0))^ - <X j (0)X fe (0))S) . (22) 

To compute the covariances in the above expression, we 
need to take a slice through the Gaussian steady-state 
distribution at the indicated fixed values of Xj. The dis¬ 
tribution of the remaining variables X' in this slice is: 


j> ss (x'|w 


ex P o E G Ik X :X k 


j,k^i 




= exp - 9 E^E^-E^E 


(23) 


(24) 


jk 


The mean of the new distribution depends on the vec¬ 
tor fj,j = CT 1 X, and the new covariance matrix C' = 
<X,-(0)X fc (0)>^ is found by removing the zth row and 
column from C -1 and then inverting it. The key prop¬ 
erty of this distribution for our purposes is that C' is 
independent of the value of Xj. It does depend in gen¬ 
eral on our choice of which row and column to remove, 
but it does not change when we vary Xj from 0 to some 
other value. Applying this fact to equation (22), we see 
that the term in parentheses equals zero for all values of 
Xj, and so equation © becomes: 


provided can extend the range of validity for this repre¬ 
sentation beyond the linear-response regime to which it 
has been traditionally restricted. 

This equation also looks similar to the result of macro¬ 
scopic fluctuation theory that demonstrates the equality 
of the log of the steady-state distribution for the den¬ 
sities and the “excess work” El- There is no obvious 
relationship with this result, however, because we have 
defined A(W)(Xj) by subtracting a constant from the 
bare mean work, rather than decomposing the current 
into symmetric and asymmetric parts. 


B. Departure from Linear-Response J aa (E) 

In the one-dimensional case, where X = J — J ss , the 
dependence of J SS (E) on E is fully determined by the mi¬ 
croscopic reversibility assumption underlying (251 com¬ 


bined with the linear equation of motion (16) and the 
definition of work (18). The result is 


J «{E) = ^<J 2 ) eq . (26) 

We can compare this to the prediction of the Green- 
Kubo formula of near-equilibrium linear response theory 
(cf. i): 

/•OO 

J^(E) = pVE / dt(J(t)J(0)) eq . (27) 

Jo 

Under the linear overdamped Langevin dynamics ana¬ 
lyzed above, this becomes: 

4 R (£) = W) {j2)e * (28) 

we see that the system departs from the linear-response 
regime when the damping rate A begins to depart from 
its equilibrium value A(0). The size of the departure is 
given by 


J ss - 4 R _ -A(O) 

ME) ■ 


(29) 


Equations (251, (|T6|) and (18) also impose a definite 


Inp ss (X'j) = —/3[F(Xj) — A(W)(Xj)] + A f. (25) relationship between A and B: 


For this strictly linear case, then, the A {W n ) c vanish 
for all n > 1. This remains true even if the dynamics 
are not Markovian for the current J on its own (as is 
clearly the case in panel (a) of Figure [2] below), as long 
as there exists a set of observables X with sufficiently 
many additional degrees of freedom that dynamics of the 
form given in equation ( fToj ) apply. 

Equation (25) is equivalent to the McLennan ensemble 


B(E) 2 
2 A(E) 


= M 2 )e q . 


(30) 


When E = 0, this is an analogue to the Einstein relation, 
constraining the relationship between the effective “mo¬ 
bility” and effective “diffusion coefficient” for the fluctu¬ 


ations of the current in equilibrium. Equation (30) says 


as presented in equation (3.13) of [8] (our averages are 
evaluated under the driven dynamics, instead of the un¬ 
driven dynamics, but that change is 0(e 2 ) in their nota¬ 
tion) . But as we will show below, the derivation we have 


that this relationship continues to hold beyond equilib¬ 
rium, even where equation (29) indicates a large differ¬ 


ence between the actual current and the predictions of 
linear response theory, as long as the equation of mo¬ 
tion remains linear. This relationship guarantees that the 
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variance of the steady state distribution remains equal to 
its equilibrium value in this regime, even as the external 
drive alters the rate of relaxation. A related constraint 
also exists in the multi-dimensional case, where more pa¬ 
rameters come into play, and the whole covariance matrix 
of the steady-state distribution must remain equal to the 
equilibrium matrix. 


C. First Correction 


We now allow nonlinear terms to be introduced into 
the coarse-grained equations of motion. Our first goal is 
to obtain an expression in terms of work and temperature 
for the dimensionless expansion parameter that measures 
the significance of the nonlinearities in the computation 
of the steady-state distribution via equation (131. We 


also seek to determine which terms from the cumulant 
expansion are present at first order in this parameter. For 
this calculation, we again focus on the 1-D case, where 
the dynamics of the current J are Markovian. We begin 
by adding a small nonlinear term added to the linear 
dynamics for X = J — J SS (E): 


X = A(E)X + h -X 2 + B(E)£(t) 


(31) 


where £(t) is again a Gaussian white noise term with 
mean 0 and autocorrelation function (£(0)£(t)) = S(t). 
Note that the sign on A is positive, because we need 
to average over trajectories whose final conditions are 
specified, as opposed to the usual initial conditions. The 
coefficient /j has dimensions of 1 / [time] [current], and part 
of our task is to convert it into a dimensionless expansion 
parameter. 

We can write the solution as a power series in /x: 

X(t) = I< 0 > (t) + fiX^ (t) + n 2 X^ (*) + .... (32) 


Plugging this in to equation of motion and collecting 
terms in powers of fi, we find 

X (0) = aX (0) + B£(t) (33) 

i (1) =al (1) + ^(l (0) ) 2 . (34) 

The first is the same as the equation we solved above, 
giving 

* (0) (i) =e M (X(0) + f(t)) (35) 


with /(f) = — f t ° e~ At Bt;(t')dt'. We can use this result 
in the second equation and solve in a similar way to ob¬ 
tain 

X (1) (t) = - [° dt ' ( X(0) ^)) 2 . (36) 

Jt 2 

Integrating the perturbative solution X(t) = X^°\t) + 
//X (1) (t) +0(M 2 ), we obtain the work done for a given 


realization of the noise: 
/■O 


W=VE 


dt 


I-At. 


e Ai X(0) + e M f(t) 


VE 

~aT 

r 0 


+ 0(fx 2 ) + J ss (37) 

x (°) + ^(O) 2 + VE dt{e At f(i) 

dt'e A ^(iiX(0)f(t') + |/(t') 2 ] + 0(M 2 ) + J ss } 

(38) 


Now we need to use this expression to compute the 
A (W n ) c (X). Since /(f) is independent of the choice of 
ending state X(0), we find 


A(W)(X) 


VE v VE 
~A(E) X + A(E) 2 


X 2 + 0(n 2 ) 


A (W)W(X) + + 0(v?), 

(39) 


where A(W)^°^(X) is computed under the linear dynam¬ 
ics alone. The quadratic term in equation (39) pro¬ 


vides a correction to the equilibrium variance at order 
/x, which also breaks the relationship (30) between A(E) 
and B(E). 

To compute the higher cumulants, we first note that 
the first two terms in equation (38) are not random vari¬ 


ables, and serve as a constant offset that makes no con¬ 
tribution to the cumulants of order 2 and higher. Fo¬ 
cusing on the term in curly brackets, then, we see that 
the fi/(f') 2 term can only be part of an A'(0)-dependent 
term in A (W") 0 when it is multiplied by some nonzero 
power of the fiX(0)f(t') term. It therefore contributes 
only to the 0(ft 2 ) part of the expression and to the over¬ 
all normalization. The remaining part of the work can 
be expressed as a sum of independent Gaussian random 
variables £(f), so it is itself a Gaussian random variable. 
This implies that it has no nonzero cumulants beyond the 
variance, so that all higher cumulants are 0(fj, 2 ). The 
variance is given by 


(W 2 )t^ x = (W 2 ) ss ^x - (W)l^ x 

= 2 V Z E^ I dt 

J —oo J —» 

x(mf(t"))+o^ 2 )+Af 


2^2.. I / dt , / dt „ e A{t+t’+t") X 

— OO j t' 


(40) 


where we have absorbed the A-independent terms into 

AT. 

To simplify this, we must compute the autocorrelation 













function of f(t). We first consider the case \t"\ > |t|: 


{fit) fit")) = B 2 ds J° du e ~ A <*+“) (mfiu)) 

r o ft 

' ds due~ A(s+u) 5is-u) 

Jt Jt" 


= B 
+ B 2 
= B 2 [ 


/ 0 pO 

ds J due^^+^Sis — it) 

-1). (41) 


ds 


If \t"\ < |£|, this becomes 


A: 


B 2 


d2 

2At 

2A { 


(fit)fit")) = ^~ 2At " !)■ 

Combining the two answers, we find 

{fit)fit")) = ^ A ie- 2M ™- 1 ) 


(42) 


(43) 


where t m is equal to whichever of t , t" has the smaller 
absolute value. 

Plugging this in to equation (401, we find 
B 2 


(W 2 ) c x =2 V 2 E 2 nX 
r o 


2 A 
o 


/ U rU f U 

dt / dt’ / dt" ( e ^(*'H*-*"l) 

-oo J — oo J t' 


— e A ( t+t ' +t ")) + 0(p 2 ) +Af 


=W 


v 2 e 2 xb 2 

A 4 


0(/i 2 )+AC 


(44) 


Thus we obtain the next term in the cumulant expan¬ 
sion in equation (13): 

^A(W 2 Y{X) = ^{{W 2 Y SS ^ X (W 2 )^ 0 ) 


P 2 v 2 e 2 b 2 v , 2 


= m 


2A 4 


X + 0(/z 2 ). (45) 


By the argument given above equation (40), the rest of 
the terms in the cumulant expansion are guaranteed to 
be 0(/i 2 ). To obtain the correct dimensionless form of 
/i, we need to compare this term to the fJA{W) term, 
and see what combination of parameters controls their 
relative size. Using equation (39), we find 

yA {W 2 YiX) = p^^/3A(W)iX) + 0( A t 2 ), (46) 

so that ft = v *™* is the appropriate dimensionless 
version of ft that controls how quickly the expansion con¬ 
verges. 

We can interpret ft thermodynamically by comparing 
to the quadratic 0{p) term in (39) that corrects the 
variance of the distribution. Using the fact that the 


variance in the unperturbed steady-state distribution is 
a\ = B 2 /2A (cf. eq. 3.8.74 and 4.3.23 in jl9j), we can 
write 


ft = 4/3[A(W)(a x ) - A(W)(°V*)] (47) 


which is four times the typical extra mean work difference 
due to the nonlinear term in the dynamics, divided by 
k B T. 

Thus we conclude that the first terms in the expansion 
of p ss iX) about the linearized dynamics in this 1-D case 
are given by 


In p S siX) = -p[F{X) - A(W)(A)] - L-A(W 2 ) C (X) 


+ 0{ft 2 ) +Af 


(48) 


with f, as defined above. 

In the multidimensional case, defining the ft explicitly 
in terms of the model parameters is more challenging, 
but it is reasonable to suppose that the thermodynamic 
expression (47) should still give a good estimate of its 


size. In the next section, we will apply this result to 
our numerical simulation of a sheared colloid, which will 
require the finite-time approximate form 


lnp ss (X) » -P[F(X) - <W) At pO] 
+ 0{ft 2 ) +J\T 


|V 2 >A t m 

(49) 


which holds when A(> t = 1 /A, with the finite-time 
cumulants defined in equation (14). 


IV. SHEAR THINNING EXAMPLE 


Strongly driven colloids are commonly used as exam¬ 
ples of far-from-equilibrium steady state systems where 
the relationship between currents and applied fields can 
violate the predictions of linear response theory (cf. [20b 
[24]'). Using the shear stress a xy as the current and minus 
the shear rate —7 as the conjugate field, one can describe 
departures from linear response in terms of the depen¬ 
dence of the viscosity 77 = —< 7 X y/'i on 7. In the linear re¬ 
sponse regime, a xy is proportional to 7, and 77 is constant. 
As the shear rate is increased in a typical colloid, a xy (7) 
becomes sublinear, so 77 decreases and the suspension is 
said to “shear thin” (cf. [23]). In this section, we describe 
how to measure (W)At{&xy) and ^f(W 2 ) Xt {a xy ) in com¬ 
puter model of a colloid that exhibits shear thinning. We 
will show that the above expansion in the degree of non¬ 
linearity converges quickly in this case, so that the ft —>• 0 
form given in equation ( [25| generates a qualitatively cor¬ 
rect description of the shear thinning phenomenon, while 
the 0(/t) correction in equation ( |49| is sufficient to main¬ 
tain quantitative agreement with the actual distribution 
well into the thinning regime. 
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FIG. 1. (Color online) Shear cell with periodic boundary con¬ 
ditions along flow direction. The reflecting walls on the top 
and bottom are separated by a distance d. The top wall moves 
at constant speed v in the x direction, while the bottom wall 
is fixed, causing a linear gradient 7 = v/d in the solvent flow 
velocity along the y direction. 


A. Setup 

Consider a suspension of small identical spheres in a 
Newtonian solvent at a fixed temperature. The particles 
are small enough that Brownian motion can equilibrate 
their spatial configuration rapidly compared with the du¬ 
ration of a typical experiment, producing a steady state 
independent of initial conditions. Electrostatic repulsion 
keeps the spheres far enough apart that the disturbance 
each particle creates in the flow field has no effect on the 
trajectories of the other particles, while ions in the sol¬ 
vent screen the charges and exponentially suppress the 
interaction at large separations. 

A nonequilibrium steady state can be created by mov¬ 
ing one wall of the chamber containing the suspension 
at a constant velocity v while keeping the opposite wall 
fixed, thus setting up a steady shear flow in the gap of 
width d between the walls. (The wall velocity is ulti¬ 
mately due to some time-varying fields A (t) - like the 
fields inside an electric motor). The strength of the shear 
flow can be quantified in a form independent of the sys¬ 
tem dimensions as the “shear rate” 7 = v/d. A constant 
shear rate can be maintained by using periodic boundary 
conditions in the flow direction (which can be approxi¬ 
mated in experiment by using a cylindrical geometry). 
As indicated in Figure [I] we will define coordinates such 
that the moving wall travels in the +x direction and the 
y axis points from the stationary wall to the moving wall. 

Two important dimensionless parameters for the dy¬ 
namics of a sheared colloid are the Reynolds number Re 
= p^a 2 /po and the Peclet number Pe = 7 a 2 b/ksT. Here 
p is the mass density of the fluid (assumed to be compa¬ 
rable to the density of the particles), a is the radius of a 
particle, pg is the viscosity of the suspending fluid, b is 
the drag coefficient of a particle (= 6npoa for a sphere 


with no-slip boundary conditions), ks is Boltzmann’s 
constant, and T is the temperature of the heat bath cou¬ 
pled to the fluid. Re measures the importance of inertia 
relative to viscous drag, and Pe measures the importance 
of motion by convection in the shear flow relative to diffu¬ 
sive motion in a dilute suspension (when the suspension 
becomes sufficiently dense, the equilibrium relaxation is 
slowed down by the interactions among the particles, 
and the relevant dimensionless parameter becomes the 
density-dependent Weissenberg number). Pe thus mea¬ 
sures the distance from equilibrium, so that Pe 1 gives 
rise to the linear-response regime, and Pe 1 constitutes 
the “far from equilibrium” regime where shear thinning 
occurs. 

In the overdamped limit where Re <C 1, the instanta¬ 
neous velocity of the particles can be regarded as fully de¬ 
termined by their spatial configuration (up to the rapidly 
equilibrating contribution from Brownian motion), so the 
set of particle positions is sufficient to define the full mi¬ 
crostate, and the corresponding simulation algorithm be¬ 
comes very simple. Re can be kept in this regime while 
sweeping Pe up to any desired maxim um value Pe max by 
choosing a viscosity such that po \J pPe max fcBT/a. 

B. Equations of Motion 

To describe this system mathematically, we will use 
the model employed in [23J [23] for the investigation of 
departures from near-equilibrium linear-response behav¬ 
ior in nonequilibrium steady states. This model can be 
numerically simulated with the dilute limit of the Brow¬ 
nian Dynamics of Ermak and McCammon |26| or of the 
Stokesian Dynamics of Brady and Bossis where hy¬ 
drodynamic interactions are ignored, and the Re -C 1 
limit is invoked so that the particle inertia becomes ir¬ 
relevant. These assumptions lead to the following dis¬ 
cretized equations of motion for the position ( Xi,yt ) of 
particle i confined to move in two dimensions: 

Xi(t + At) = Xi{t) + yi(t)pAt + ^ ^ x • FjiAt + Ax r t 

(50) 

Vi{t + At) = yi (t) + * y ■ F :i At + A v r i ( 51 ) 
3 

where b is the drag coefficient for the particles, and F ? , 
is the force exerted on particle i by particle j. We choose 
the force to be a screened Coulomb repulsion, with poten¬ 
tial energy U(r) = kBTe~ rA zlB/r as a function of the 
distance r separating a pair of particles. A is the screen¬ 
ing length, Ib is the Bjerrum length, and z is the number 
of elementary charges on each particle. Ax\ and A y[ are 
random displacements due to Brownian motion. Brow¬ 
nian and Stokesian dynamics assume that the solvent 
degrees of freedom always equilibrate before the spatial 
configuration of the colloidal particles can change appre¬ 
ciably, so that the Brownian displacements can be chosen 
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from a Gaussian distribution whose variance (2 /b)At 
is related to the temperature and drag coefficient by the 
Einstein relation. Equations ([501) and (51) are then sim¬ 


ply iterated with a small enough time step that the re¬ 
sults are insensitive to variations in time-step size }28j . 

As mentioned at the beginning of this section, we 
are considering the case where the particle size is much 
smaller than A or zIb , so that “hydrodynamic interac¬ 
tions” (particle-particle interactions mediated by distur¬ 
bances in the solvent flow) make a negligible impact on 
the particle trajectories. This is what allows us to use the 
“dilute limit” of the Stokesian or Brownian Dynamics, 
where the mobility and resistance tensors are diagonal 
and independent of particle positions. Another conse¬ 
quence of this limit is that the actual particle radius a 
does not appear in the equations of motion; we therefore 
use the screening length A as the microscopic length scale 
for computing the Peclet number and measuring distance 
from equilibrium. 


shear rate 7 . The remaining term a xy is independent of 
the particle positions, so the work done by the moving 
wall will depend on the particle configuration through 
aL,. alone. 

When the shear rate is small compared to the diffusive 
relaxation rate, the overall viscosity of the suspension can 
be computed from the equilibrium fluctuations in <j xy us¬ 
ing linear response theory [3j . As the shear rate continues 
to increase, the viscosity begins to deviate from this value 
as the suspension shear thins. In the following sections, 
we will use equation (49) to determine the most probable 


value of <j xy and hence the contribution 07^/7 to the vis¬ 
cosity in both the linear response regime and the shear 
thinning regime. We start in section |TVD| by demonstrat¬ 
ing how to extract (W)a t(cr xy ) and ^(W 2 ) c At (a xy ) from 
the simulation in a way that should also be experimen¬ 
tally accessible. Then in section [V] we plot the results 
over a range of shear rates, and compare the prediction 
of equation (49) with the observed steady-state distribu¬ 
tion. 


C. Shear Stress 


The macroscopic viscosity of the whole suspension at 
equilibrium will be larger than 770 , because both the dis¬ 
turbance of the flow field produced by individual particles 
and the mutual repulsion between pairs of particles make 
the suspension harder to shear than the bare fluid. As 
the suspension is sheared, however, the contribution of 
the particle repulsion to the viscosity decreases, and the 
suspension shear thins (cf. (25]). The particles cause the 
shear stress to vary with position in the suspension, so 
we define an overall shear stress for the system by av¬ 
eraging the local shear stress at the moving wall of the 
system over the whole wall area. This will be convenient 
for computing the work done by the moving wall later 
on, and gives us a macroscopic parameter that can be di¬ 
rectly observed in experiment via a measurement of the 
force applied to the wall. As shown in the Appendix, 
for a suspension of particles in a Newtonian solvent in 
the limit of zero Re with no hydrodynamic interactions, 
the instantaneous mean shear stress a xy u exerted by the 
fluid on the moving wall is: 

<y U = °ly + °xy ( 52 ) 

The first term depends on the force F ly exerted by each 
particle i on the other particles j: 

17 *y = • (53) 

i¥=J 

Here V is the system volume, x is the unit vector in 
the +x direction, and A— y,. The right-hand 
side can be unambiguously determined from the system 
microstate, which we are taking to be the list of positions 
of all the particles. We can therefore choose X = a xy as 
our coarse-grained observable and test how well the first 
terms in our expansion describe its fluctuations at a fixed 


D. Measuring the Work Cumulants 


The rate at which the moving wall does work on the 
fluid is just the force —Aa xy n it exerts against the fluid 
(where A is the surface area of the wall) times the speed 
of the wall yd. Using equation (]52|), we thus obtain: 


W = -Vi aiy + W 0 . 


(54) 


where Wo is the part of the work that does not depend on 
the configuration of the particles. Since Wo only affects 
the normalization, but not the shape of the distribution, 
we will set it to zero for the purpose of the calculations 
in this section. Then we can treat erL, as the current J 

y 

and minus the shear rate —7 as its conjugate field E. 


Using equation (54), we can now compute the work 
done along any stochastic trajectory <J xy (t) by simply in¬ 
tegrating the trajectory with respect to time. We can 
then compute the distribution of W for trajectories end¬ 
ing at a given a xy value by letting the system relax to the 
steady state at some value of 7 and run there for a long 
time, while continuously recording the fluctuations in a xy 
(which can in principle be determined directly from the 
fluctuations in the force applied to the moving plate in 
an experiment). As shown in Figure ([ 2 ]), we then choose 
some time interval At longer than the timescale r of re¬ 
laxation to the steady state, and compute both the work 
and the final value of a xy for every segment of length At 
in the whole trajectory. Finally, we bin the work out¬ 
puts by the corresponding final value of a xy to obtain 
the distribution of work for each bin, from which we can 
compute the cumulants (W" )a 

Figure [ 3 ] shows the zeroth order term (W)At(&i y ) 


in the expansion (491 and the first-order correction 


03/2 )(W 2 ) A t( a xy) as a function of ending value of oy 
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FIG. 2. (Color online) (a) cr xy (t) averaged over all trajectory 
segments of length At = 0.5 that end at a specified value, 
from a simulation run with the same parameters as Figure [4] 
and Pe = 19. Two arbitrarily chosen ending state restrictions 
are indicated by dotted lines, with the corresponding average 
trajectories shown in the same colors, (b) Measured proba¬ 
bility distributions of work done over all the trajectories that 
go into the averages of panel (a), shaded in the same colors, 
(c) A portion of the raw a xy timeseries from which the other 
two panels were generated. The ending state restrictions of 
panel (a) are indicated with dotted lines of the same color, 
and two typical trajectory segments ending at these values 
are shaded in these colors. The shaded areas are proportional 
to the interaction-dependent part of the work, according to 
equation (54). 


at three different values of the dimensionless shear rate 
Pe. 


V. SIMULATION RESULTS 


With this procedure in hand for extracting (W)a* and 
(W 2 ) < kt( a xy)i we can use equation (49) to compute the 
steady-state distribution of a T xy from the work statistics, 
and compare it to the directly measured distribution in 
our simulation. 

We simulated a sheared colloidal monolayer of N = 100 
particles using the equations of motion (50) and (51 1 [5B]. 
The colloid was confined to a square box of side length 20, 
with reflecting boundary conditions on the moving wall 
and the opposite wall, and periodic boundary conditions 
on the other sides. This concentration is sufficiently di¬ 
lute that the equilibrium relaxation time does not vary 
much with changes in concentration, so the Peclet num¬ 
ber should still be the relevant parameter for measuring 
distance from equilibrium. The other parameters were 
chosen as ksT = b = A = zIb = 1. We ran this simula¬ 
tion for 20 different values of Pe, from 0 to 19, generating 
trajectories with lengths up to t = 120, 000 in the given 
units, with time step size 0.001. The simulations were 
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FIG. 3. (Color online) Using the method illustrated in Figure 
[ 2 } we compute (W)At(al y ) and {P/2)(W 2 ) c At (a I xy ) for a range 
of values of a xy , and plot this data for Pe values 4, 10, and 
19, increasing from bottom to top. Also plotted is the free 
energy F extracted from the Pe = 0 simulation run. 


initialized with uniform random distributions of particle 
positions, and the initial transients were removed from 
the timeseries before analysis. 

To apply equation (49), we first need to find the equi¬ 
librium free energy F(a^ y ) as a function of a T xy . This can 
be extracted from the simulation by extracting the dis¬ 
tribution of equilibrium fluctuations from a run at 7 = 0 , 
taking the natural logarithm, and multiplying by —ksT. 

Figure [3] shows how F and the other two terms in equa¬ 
tion ( |49| ) depend on <J xy , with the nonequilibrium terms 
evaluated at three different values of the shear rate. F 
is parabolic near a xy = 0 , but requires a fourth-order 
polynomial to fit the far tails. {W )At starts out linear 
at low shear rates, starts curving slightly by Pe = 10, 
and becomes noticeably quadratic by Pe = 19, indicat¬ 
ing that the 0(jl ) term has become important. (W 2 ) At 
is independent of a xy at low shear rates, but starts be¬ 
coming <j^ y -dependent at about the same shear rate as 
(W)a t begins to deviate from linearity, as expected from 
our analysis in section |III C| 

Figure [4] shows the location of the peak of the steady- 
state distribution as a function of Pe computed using 
equation ( |49| ), compared with the distribution directly 
sampled from the simulation. We have also plotted the 
result ignoring the A {W 2 ) c term, to show the size of the 
impact of the O(p) correction compared to the zeroth 
order expression (25). When we compare both curves 
to the most probable values of a xy actually measured in 
the simulations, we see that the zeroth-order expression 
correctly captures the qualitative shear thinning behav¬ 
ior: a xy departs from its linear-response dependence on 7 
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FIG. 4. (Color online) Top: Location aiy of peak of probabil¬ 
ity distribution, computed via equation (491, compared with 
directly measured most frequent values in the simulation [2S]. 
Also plotted is the prediction using the mean work alone via 
equation (251, without the O(jl) variance term, as well as the 
linear response prediction Bottom: Relative size of 

departure from linear-response prediction. 


and eventually saturates, causing the interparticle force 
contribution —cr^ /7 to the viscosity to fall off as 1 / 7 . 
This approximation predicts that the saturation occurs 
sooner than it actually does, but the first-order term ap¬ 
pears to entirely compensate for the discrepancy. The 
straight line is the linear-response prediction for the mean 
shear stress, computed from the equilibrium fluctuations 
using the Green-Kubo formula given in equation (27). 
The bottom panel shows the relative difference between 
the linear-response prediction and the location of the ob¬ 
served peak of the prob abili ty distribution, as defined in 
equation (l29|) of section III 


In Figure [5] we have plotted the normalized probabil¬ 
ity distributions at Pe = 4, 10, and 19, showing the di¬ 
rect measurement and the prediction of equation (49), 
smoothed with polynomial fits to the data in Figure [3J 
Also plotted is the equilibrium distribution shifted to the 
location of the new peak, to better visualize the change 
in the shape of the distribution at high shear rates. The 
variance is still within 2 % of the equilibrium value at 
Pe = 4, even though the relative deviation of the mean 
from the linear-response prediction has already reached 
30%. This is an example of a case where the constraint of 
equation (301 derived from the linear approximation re¬ 
mains in effect beyond the linear-response regime. Equa¬ 
tion (491 successfully accounts for the change in variance 
visible at Pe = 10, although it fails to fully capture the 




FIG. 5. (Color online) Full probability distribution from equa¬ 
tion ( |49| ), compared with shifted equilibrium distribution and 
with distribution directly sampled from simulation at Pe = 4, 
Pe = 10 and Pe = 19. 


asymmetry that enters the distribution at Pe = 19, which 
should depend on higher-order terms according to the 
analysis of section (III C). 


VI. DISCUSSION 


We have shown that the probability distribution of a 
macroscopic observable in a driven steady state can be 
written as a perturbation expansion in the nonlineari¬ 
ties of the coarse-grained dynamics, which can converge 
quickly in some systems even under strong driving. The 
approximate formulas (25) and (48) give a good descrip- 
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tion of the statistics of the steady state under the follow¬ 
ing conditions: 

• The transition rates between system states must 
satisfy the “microscopic reversibility” condition 
given in equation ([I]). 

• The system exponentially loses memory of its ini¬ 
tial conditions, with a decay time that is short com¬ 
pared to the duration of the experiment or simula¬ 
tion. 


• The contribution to the mean work difference 
A (IT) (defined in equation ©) due to nonlinear 
terms in the coarse-grained fluctuation dynamics is 
small compared to the thermal energy scale ksT. 


The simulation results presented in Section IV D and 
Section [V] illustrate the physical meaning of the terms in 
the expansion, and confirm that it can converge rapidly 
even when external drives are strong enough to make the 
peak of the steady-state distribution depart significantly 
from the linear-response prediction. 

Evaluating the distribution by measuring (IT)At and 
(W 2 )! t numerically, as done in section [v] is much more 
computationally expensive than a direct measurement of 
the distribution from the simulation data. This compu¬ 
tation was done in order to demonstrate that the statis¬ 
tics of the simulated system can be captured by the first 
terms in the perturbation expansion, but the real signif¬ 
icance of our result lies in the physical intuition that can 
be obtained for the steady-state behavior of systems that 
reside in this regime. 

The zeroth-order expression (251 constitutes a natural 
extension of the first attempts by Einstein and his con¬ 
temporaries at understanding macroscopic fluctuations 
(cf. [29]). They interpreted the equilibrium macroscopic 
distribution to mean that the probability of a fluctuation 
in a given observable is exponentially suppressed in the 
ratio of the work that must be done by the rest of the de¬ 
grees of freedom in the system to produce this fluctuation 
to the thermal energy scale fc^T. Equation (25) simply 
incorporates the fact that in a driven system, some of 
the work can be done by the external drive instead of by 
other internal degrees of freedom. By subtracting off the 
extra work done by the drive on the way to a given fluc¬ 
tuation (compared to the work done on the way to a fixed 
reference state), we account for the “help” provided by 
the drive in supplying the work required to reach each 
state. The first correction term from the perturbation 
analysis begins to account for the stochasticity in this 
extra work: even if two possible fluctuations extract the 
same amount of work on average from the drive, one can 
be less likely than the other if its distribution of extracted 
work extends further towards zero, allowing the fluctua¬ 
tion to be reached by paths that receive little help from 
the drive. 

For strictly linear relaxation of a single dissipative cur¬ 
rent, the mean extra work from the drive can be com¬ 
puted up to an additive constant C in terms of the re¬ 


laxation time t(E) = 1/A(E) using equation (20) from 


section III A(W)(J) = VEJt(E)+C. For systems close 
to this regime, the nonlinear response of the steady-state 
distribution to an increase in the drive can thus be esti¬ 
mated from knowledge of the behavior of the relaxation 
time. In our sheared colloid case, we can understand the 
shape of the stress vs. shear rate curve in Figure [4] this 
way. The equilibrium relaxation time r(0) is set by the 
diffusive time scale 6 a 2 /fcsT, but in the 7 —» 00 limit, 
t(E) should be dominated by convective “stirring” from 
the shear flow with time scale I/ 7 . The shift in the peak 
of the cr^y distribution should thus increase linearly with 
the strength of the conjugate field E = —7 near equi¬ 
librium, with the expected linear-response coefficient via 
equation (28); but the peak should stop shifting once 


t( 7 ) reaches its asymptotic I /7 behavior, which cancels 
out all the 7 -dependence. 


As we pointed out in equation (30) in section |III| the 
linear regime also exhibits the remarkable property that 
the relaxation time r(E) and the coefficient B(E ) on the 
noise term in the Langevin equation are related by a gen¬ 
eralized Einstein relation or fluctuation-dissipation theo¬ 
rem. A decrease in r( 7 ) must accompanied by an increase 
in the strength of the noise term in the coarse-grained 
dynamics, in such a way that the variance of the steady- 
state distribution remains exactly equal to the variance 
of the equilibrium distribution. The size of the fluctua¬ 
tions thus remains unchanged in the regime of linear re¬ 
laxation, even as the dynamics are significantly altered, 
with a new time scale and (as in the case of this sheared 
colloid) a new underlying physical mechanism. The first 
panel of Figure [5] confirms that this constraint remains 
applicable even when the linear-response prediction for 
the mean shear stress is no longer valid. 


One could now explore the application of the theory 
to richer phenomena, such as hydrocluster formation in 
shear thickening colloids (cf. [5(]J 31]), where the distri¬ 
bution of some measure of typical cluster size could be 
computed in terms of the mean work done for trajec¬ 
tories that end at a given value of that parameter. It 
would also be interesting to investigate whether suspen¬ 
sions of active particles in a quiescent solution can be 
analyzed in this way, possibly delivering further insight 
into the “freezing by heating” transition that takes place 
at a critical propulsion rate in both experiment and sim¬ 
ulation [3211341 . 


Equation (13) should also be readily generalizable to 
chemical as opposed to mechanical driving, using the ex¬ 
tensions of equatio n ([l] ) to chemical reaction networks 
mentioned in section pljjl I HT(7| . Once the quantitative re¬ 
lationship between the bulk quantities of interest and the 
work rate are understood, this generalized result could 
shed light on the steady-state properties of biologically 
relevant systems, such as active actin-myosin networks 
driven by ATP hydrolysis 
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Appendix: Derivation of Mean Wall Stress 
(Equations |52| and 53) 


Formulas for determining the particle contribution to 
the shear stress of a colloidal suspension have been known 
for a long time, and received an especially careful treat¬ 
ment by G.K. Batchelor in the 1970’s [3UI3H1. The estab¬ 
lished literature mainly deals with the mean shear stress, 
either averaged over an infinite ensemble of systems or 
over an infinitely large system. The statistical unifor¬ 
mity of the system can then be invoked to argue that 
the mean stress over a typical 2-D slice through the sys¬ 
tem is equal to the mean stress averaged over the whole 
system volume. Although the wall is not a typical 2-D 
slice, because the boundary condition modifies the parti¬ 
cle distribution, the fact that there is no mean net force 
on any part of the system when it is in steady state im¬ 
plies that the mean stress on all parallel 2-D slices must 
be the same. The average over an infinite system vol¬ 
ume must therefore also be equal to the average over an 
infinite wall [38i |. 

For the purpose of this paper, it is not enough to know 
the ensemble- or infinite-system-averaged mean. We need 
to look at the fluctuations about the mean in order to ap¬ 
ply our procedure for empirically determining the mean 
renormalized work and the equilibrium free energy as a 
function of the shear stress. Therefore we need to go 
back through the derivation, and examine the instanta¬ 
neous value of the shear stress at the wall in a suspension 
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of a finite number of particles. 

In this appendix, we prove that the instantaneous 
shear stress exerted by the fluid on the moving wall of 
the shear apparatus described in section III, averaged 
over the moving wall area, is 

_wall _ J | _0 /a i\ 

® xy ®xy ®xy 

where a 1 is defined by 

°iv = ( A - 2 ) 

and <r xy is independent of the particle positions. 

We start by giving some necessary background on the 
behavior of shear stress in low-Re Newtonian fluids. To 
make this proof accessible to readers less familiar with 
hydrodynamics, we then map to a mathematically anal¬ 
ogous problem in electrostatics (which turns out to be 
a homework problem from Griffiths Electricity and Mag¬ 
netism 00 problem 3.44a]). After presenting the solu¬ 
tion to this electrostatics problem, we finally map back 
to hydrodynamics to obtain our final result. 


1. Stress in Newtonian Fluids 

The shear stress a xy is an off-diagonal component of 
the 3-by-3 stress tensor a. a is defined at each point in 
the fluid such that n • a is the force per unit area exerted 
from below on a surface element at that location with 
unit normal vector n. By “from below,” we mean from 
the side opposite to the direction of the normal vector. 
We will focus on the x column cr ■ x to obtain a vectorial 
quantity that will be easier to visualize. 

By the definition of the stress tensor above, the x- 
component of the force on a region Cl of fluid is given 
by 


F, r = - 


dA ■ a ■ x 


id n 


= — dW ■ a ■ x 
J n 


(A.3) 

(A.4) 


which implies 


cr • x = -r] 0 Vu x (A.6) 

where u x is the cc-component of the fluid velocity field, 
and 770 is the (constant) viscosity of the solvent. Com¬ 
bining this with the previous equation gives us the set of 
equations 

?7o V 2 u x = f x (A.7) 

a • x = -77 0 Vii x (A.8) 

that together fully determine cr-x for a given set of bound¬ 
ary conditions. 


2. Mapping to Electrostatics 


Equations (A.7) and (A.8 ) suggest a mapping to elec¬ 


trostatics. T]ou x is the analog to the electric potential </>, 
a ■ x is the analog to the electric field E, and —f x is the 
analogue to the charge density p. With these mappings, 
the mathematics of the problem are identical to electro¬ 
statics, and we can do everything in terms of E, <p and p 
until we map back at the end. 

The only remaining piece of setup is to map the bound¬ 
ary conditions and the “charge distribution.” The non¬ 
slip boundary condition requires that every part of the 
fluid in contact with a non-rotating rigid surface must 
share the same velocity. Since the electric potential <f> 
is the analog of the ^-component of velocity, this implies 
that non-rotating surfaces behave like perfect conductors 
- they are always equipotentials. In particular, the con¬ 
straint that the bottom wall is fixed and the top wall 
moves at constant velocity v implies that the walls of the 
shear cell become parallel conducting plates separated by 
a distance d, with fixed electric potential difference A <f>. 
The problem of determining the total force on the walls 
is thus equivalent to determining the induced charge on 
these conducting plates. 

The particles, however, are allowed to rotate. Their 
boundary conditions are therefore more complicated, in¬ 
volving the other columns of the stress tensor. Specifi¬ 
cally, we have 


where dA is an infinitesimal area element of the bound¬ 
ary dCl pointing along the outward normal direction, and 
dV is an infinitesimal volume element. We add the minus 
sign because we are computing the force on this surface 
from the outside. The second line results from the diver¬ 
gence theorem. Since this holds for every possible region 
O, we conclude that the integrand of equation (A.4) is 
equal to minus the x component of force per unit volume 
f x exerted by the surrounding fluid on an infinitesimal 
volume element: 


V ■ a • x = -f x . (A.5) 

Finally, we must invoke the assumption that the solvent 
in which the particles are suspended is a Newtonian fluid, 


u = Cl x r j_ 


(A.9) 


for all points on the surface of the sphere, where rj_ is 
the vector pointing from the center of the sphere to the 
surface point, projected onto a plane perpendicular to 
the angular velocity vector Cl. Cl and the center-of-mass 
velocity u cm are free parameters that must be adj usted 
so as to be consistent with equations (A.7) and (A.8). 
The resulting restriction on a is 


a = —770 V (Slxril u cm ). 


(A.10) 


To determine the charge distribution, we use our as¬ 
sumption of low Re to require the total force on any vol¬ 
ume element to vanish. In the electrostatic analogy, this 
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implies that the solvent is uncharged, and all charge must 
reside at the walls or on the particles. The interparticle 
repulsion exerts force on each particle that must be can¬ 
celed by the friction of the fluid in order to satisfy the 
requirement of zero total force. This implies that the to¬ 
tal “charge” on each particle must be qi = YYj^i F ji ' x, 
where F :]1 is the force exerted on particle i by particle j. 
The distribution of this total charge over the surface of 
each sphere is not fixed in advance, however, and must be 
determined by solving equations (A.71 and (A.8) (along 
with the corresponding equations for the other compo¬ 
nents of the stress tensor) with the boundary conditions 
just described. The decision to “ignore hydrodynamic 
interactions” mentioned in the main text allows us to 
greatly simplify the problem of determining these distri¬ 
butions, by solving the equations for each particle indi¬ 
vidually, with boundary condition E — > —(Afi/cfyy far 
from the sphere. This approximation ignores the effect 
of the other particles and of the induced wall charge on 
the charge distribution over each sphere. The solutions 
obtained under this approximation are independent of 
the particle positions, which will be important later on. 


3. Obtaining the Induced Charge on the 
Conducting Plate 


Our problem is thus reduced to determining the in¬ 
duced charge on a pair of conducting parallel plates at 
fixed electric potential due to a given charge distribution 
inside. 

We start by splitting the charge on the plates into two 
parts, following the strategy of Batchelor in his treatment 
of the effect of particle interactions on mean shear stress 
[39] . The derivation will resemble Batchelor’s in many 
ways, despite the electrostatic language, but adds the a 
new element by considering the wall stress due to a given 
instantaneous configuration of particles as opposed to an 
ensemble average of all possible configurations. 

The first part of the charge is the part required to 
maintain the electric potential difference A (f> in the ab¬ 
sence of any additional charges between the plates: Qo = 
AA(f>/d on the top and —Q o on the bottom. To find the 
remaining charge, we can solve for the case where the two 
plates are grounded. When we add up the two charge 
distributions, the resulting field is guaranteed to produce 
the desired constant electric potential difference. The 
case of grounded plates is problem 3.44a in Griffiths, as 
mentioned above, and we will follow his method to solve 
it gD]. 


tial E = — V0 to obtain 


j dVS7cj) i • E 2 = 

f dV<t> rV • E 2 

dV(j>ip2 

(A.11) 

j dVVfo • Ei = 

j" dVfoV • Ej 

dVcj)2Pi 

(A. 12) 


where we are integrating over all space, and have used 
integration by parts to switch the V from <f> to E. 

We thus obtain Green’s Reciprocity Theorem: 


d,V(j)ip2 = / dVfcpi- 


(A.13) 


Now we use this relation to compute the induced charge 
on our plates. We will start by computing the induced 
charge due to a point charge q at location r = ( x,y,z ). 
We will work in coordinates where the bottom plate is at 
y = 0 and the top is at y = d. 

To apply the Reciprocity Theorem, we choose for pi 
the actual charge distribution we are analyzing, with the 
point charge between the grounded parallel plates. We 
define Q + as the total induced charge on the top plate 
and Q- as the total induced charge on the bottom plate. 
For p 2 , we choose a charge distribution with conducting 
plates in the same locations, but with the top plate fixed 
at electric potential 4>q above the bottom one, and with 
no charge in the space between them. The LHS of the 
Reciprocity Theorem vanishes, because <pi = 0 whenever 
P 2 is nonzero. The RHS has a contribution from the 
charge distribution on the top plate, and a contribution 
from the particle. If the plates are infinite, then the po¬ 
tential a distance y above the bottom plate in scenario 
2 is exactly (y/d)cj) o- This will still be a good approx¬ 
imation in a finite system for charges that are not too 
close to the edges of the system, which will be true for 
the charges on the vast majority of the spheres when the 
number of spheres is large. Thus we obtain: 


0 — 4>oQ+ + 


(A. 14) 


Solving for Q + , we find 


Q + = ~d q - 


(A. 15) 


Griffiths starts by having the student derive a relation 
known as Green’s Reciprocity Theorem. (This theorem is 
closely related to a result due to Lorentz in hydrodynam¬ 
ics, which Batchelor employs in his analysis gT|.) Con¬ 
sider two distinct charge distributions pi(r) and p 2 (r), 
which produce electric fields Ei(r) and E 2 (r), with elec¬ 
tric potentials <^i(r) and 0 2 (r). Now use the Maxwell 
Equation V ■ E = p and the definition of electric poten- 


Now we again use the linearity of our equations to ob¬ 
tain the total induced charge by summing up the contri¬ 
butions from all the infinitesimal charge elements in the 
distribution. A convenient way to perform this sum is 
to split up the charge distribution on each sphere into 
two parts: a spatially uniform part equal to the mean 
surface charge on the sphere, and spatially varying part 
that integrates to zero over each sphere surface. 
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a. Contribution of Variations about the Mean We 
start by computing the contribution of the second part 
of the charge distribution. Since this part of the charge 
sums to zero on each sphere, every positive charge 5q has 
a corresponding negative charge — 5q somewhere else on 
the sphere. The net induced charge from each such pair 
is 


6Q + = —{y--y+) (A.16) 

where j/_ and are the coordinates of the +dq and 
— 5q charges, respectively. Now recall that by ignoring 
hydrodynamic interactions, we can solve for the charge 
distribution over each sphere without knowing its po¬ 
sition relative to the plates or the other particles. Fur¬ 
thermore, the linearity of the governing equations implies 
that the variations about the mean charge density are in¬ 
dependent of the size of the mean. This implies that the 
y-distance j/_ — y+ between any pair of charges on a sin¬ 
gle sphere is independent of the spatial configuration of 
the particles and of the total charge < 7 ,; of the particle in 
question. 

Summing over all pairs of charges from all the spheres 
in the sample, we define the quantity 

Qh = Y, 5< 3+ (A.17) 


and add them all up. We thus find 

Qi = -J2ff q *- ( A - 18 ) 

Combining the above results, we find that the total in¬ 
duced charge on the top plate is Q = Qi + Qo + Qh, with 
< 5 / the only term that depends on the particle positions. 


4. Mapping Back to Hydrodynamics 

We can now map back into the original variables (re¬ 
calling that charge is equivalent to minus the force ex¬ 
erted by the fluid) in order to obtain the total force ex¬ 
erted by the fluid on the moving wall of the shear appa¬ 
ratus: 


F wall = H 7 
i 


x • + Fq + Fh. 


(A.19) 


We can simplify this expression by using the fact that 
F = —F •• 


A W all 


1 

2d 


y x • FijAyij + F 0 + Ffj. 


(A.20) 


as the total induced charge due to the variations about 
the mean charge on the surface of the spheres. This 
quantity is independent of the particle positions, and just 
adds a constant offset to the total charge. The H sub¬ 
script stands for “hydrodynamic,” because this contri¬ 
bution comes purely from the friction of the flow field 
around each particle. 

b. Contribution of the Mean Charge To complete 
our calculation, we must compute the charge induced on 
the plate by a given configuration of uniformly charged 
spheres. Since the field of a uniformly charged sphere is 
equivalent to the field of a point charge (for points out¬ 
side the surface of the sphere), we can simply evaluate 
the point charge solution derived above for every particle, 


Finally, we can divide through by the area A of the wall 
to obtain the mean shear stress exerted on the wall by 
the fluid: 

= aly + °ly + °xy ( A - 21 ) 

where 

a xy = * ' Fjj Ayij (A.22) 

and the other two terms are independent of the particle 
positions. For notational simplicity, we combine them 
into one term in the main text, which we call a 0 ,.,,. 

> u- y 



